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ABSTRACT 

We present a new systematic way of setting up galactic gas disks based on the as- 
sumption of detailed hydrodynamic equilibrium. To do this, we need to specify the 
density distribution and the velocity field which supports the disk. We first show that 
the required circular velocity has no dependence on the height above or below the 
midplane so long as the gas pressure is a function of density only. The assumption 
of disks being very thin enables us to decouple the vertical structure from the ra- 
dial direction. Based on that, the equation of hydrostatic equilibrium together with 
the reduced Poisson equation leads to two sets of second-order non-linear differential 
equation, which are easily integrated to set-up a stable disk. We call one approach 
'density method' and the other one 'potential method'. Gas disks in detailed balance 
are especially suitable for investigating the onset of the gravitational instability. We 
revisit the question of global, axisymmetric instability using fully three-dimensional 
disk simulations. The impact of disk thickness on the disk instability and the forma- 
tion of spontaneously induced spirals is studied systematically with or without the 
presence of the stellar potential. In our models, the numerical results show that the 
threshold value for disk instability is shifted from unity to 0.69 for self-gravitating 
thick disks and to 0.75 for combined stellar and gas thick disks. The simulations also 
show that self-induced spirals occur in the correct regions and with the right numbers 
as predicted by the analytic theory. 

Key words: galaxy: disc - galaxies: evolution - galaxies: structure - methods: nu- 
merical 



1 INTRODUCTION 

The stability of gas disks plays an important role in govern- 
ing the structure of disk galaxies and in regulating their star 
formation rate. Although important insights can be obtained 
using perturbation theory (Toomre 1964, Lin & Shu 1964, 
Rafikov 2001), the onset of stability and its impact on the 
star formation and evolution of gas disks is best studied us- 
ing hydrodynamical simulations. These can follow the non- 
linear behavior of the system, which cannot be addressed 
by linear analysis. With the recent advances in computing 
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power and the development of new numerical techniques, 
we are now in a good position to treat a three-dimensional, 
isolated galaxy self-consistently. 

However, in order for a stability analysis to be mean- 
ingful and reliable, it is of paramount importance that one 
can specify equilibrium initial conditions. After all, if the 
initial disk is not in equilibrium, its relaxation during the 
first time-steps of the simulation may trigger instabilities 
that are of little relevance for our understanding of the sta- 
bility of disk galaxies. Unfortunately, no analytical solution 
is known for the density, velocity field and temperature of a 
three-dimensional gas disk in hydrostatic equilibrium in the 
external potential of a dark matter halo and/or a stellar disk. 
Consequently, previous hydrodynamical simulations have ei- 
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ther started from non-equilibrium initial conditions, or have 
resorted to iterative techniques to set-up the initial condi- 
tions, at the cost of having little control over the resulting 
equilibrium configuration. In this paper we present a new 
method that allows one to compute the density and veloc- 
ity structure of a realistic, isothermal, three-dimensional gas 
disk in hydrostatic equilibrium in an abritrary external po- 
tential. 

Hydrostatic equilibrium implies a balance between 
gravity and pressure. Gravity includes the self-gravity of 
the disk plus that of external components (i.e. dark mat- 
ter halo, bulge, stellar disk, etc), while the pressure is given 
by an equation of state p = p(p g , T), with p being the gas 
pressure, p g the gas density and T the temperature. The 
challenge is to find a p g , T and the velocity field, v, such 
that the system is self-consistent (i.e., obeys the Poisson 
equation) and in hydrostatic equilibrium. 

In the case of an isothermal, axisymmetric, perfectly 
self-gravitating disk (i.e., no external potential), the equi- 
librium disk has a sech 2 distribution (Spitzer 1942) in the 
ve rtical direction, with a scale- height that is proportional to 
sj cf I p g (R, z = 0), where c s is the sound speed. Here, cylin- 
drical coordinates, (R,cf>,z), are used to describe the den- 
sity field. This immediately shows that since p g (R, z = 0) is 
typically a decreasing function of radius, one generally ex- 
pects the scale-height to be a function of R. In particular, 
in the case of a globally isothermal disk, the sound speed 
c s oc T is constant in space, giving rise to a flaring disk, i.e., 
the scale-height increases with increasing R (Narayan & Jog 
2002, hereafter NJ02; Agertz et al. 2009). Alternatively, if 
we want to initialize a disk with a constant scale-height, a 
radial temperature gradient needs to be introduced. Tasker 
& Bryan (2006) initialize their disks to be isothermal and to 
have a constant scale-height. As indicated above, this can- 
not be an equilibrium configuration. Consequently, the disk 
is expected to experience an unavoidable relaxation process 
which makes the initialization not well-controlled and might 
potentially contaminate the physics, e.g., star formation, gas 
dynamics etc., of interest. Agertz et al. (2009) set-up their 
isothermal disks based on the local total surface density of 
gas plus dark matter. Although the scale-height of their ini- 
tial disk changes with radius, the local total surface density 
is not defined in a mathematical way and therefore elusive. 
In addition, their surface density does not follow an expo- 
nential profile. 

An important assumption underlying Spitzer's analysis 
is that the radial variation in the potential is negligible com- 
pared to that in the vertical direction. This assumption is 
supported by observation that disks typically have vertical 
scale-heights that are an order of magnitude smaller than 
their radial scale- length (van der Kruit & Searle 1981a,b). 
A well studied example is the Milky Way, whose scale-height 
is less than 200 pc for the cold gas (Jackson & Kellman 1974; 
Lockman 1984; Sanders et al. 1984; Wouterloot et al. 1990; 
see also Narayan &: Jog 2002) and roughly 300 pc for the 
stars in the Solar neighborhood (Binney & Tremaine 2008, 
Kent, Dame & Fazio 1991), compared to a radial scale- length 
of ~ 3.5 kpc. Throughout this paper we therefore follow 
Spitzer and consider disks to be 'thin', allowing us to treat 
their radial and vertical structure separately. Hence, we cau- 
tion that our method is not valid for thick disk structures. 



However, since we are mainly concerned with cold gas disks 
in this paper, this restriction is of little importance. 

Springel, Matteo & Hernquist (2005) introduce a flexi- 
ble solution for initializing a gas disk self-consistently. Basi- 
cally, they solve Eq. (2), Eq. (3) and Eq. (24) (see Section 
2) iteratively. First, they deploy a number of particles (say, 
2048 x 64 x 64) on a distorted grid structure in the radial, 
the azimuthal and the vertical directions. Unlike the live par- 
ticles, these particles are simply used as markers for mass 
distribution. Second, they compute the joint total potential 
and the resulting force field numerically with a hierarchical 
multipole expansion based on a tree code. Third, given the 
potential just computed, integrating Eq. (2) for a given mid- 
plane volume density, p g (R, 2 = 0), gives the vertical struc- 
ture of density. Fourth, adjust the midplane volume density 
to fulfill Eq. (24). Repeat the procedure between the second 
step and the fourth until the result converges. 

Although this description is quite general and flexible, 
for several reasons, this is not commonly used in the grid- 
based codes which are featured with adaptive-mesh refine- 
ment (AMR). The first and also the most fundamental one 
is that the grid structure is normally unknown before we ac- 
tually initialize the disk. Except the uniform-grid initializa- 
tion, the grid structure is automatically generated based on 
the criterion for refinement. Second, for a fully parallelized 
code, the initial data is distributed over different processors 
and memory storages. This means that the data exchange 
between processors is necessary in order to fully compute 
the joint total potential. The situation becomes even more 
technically challenging when initializing with AMR. Third, 
The vertical structure of the gas disk depends only on the 
vertical potential difference (see Eq. (7) and Eq. (9) below). 
A description of the equatorial potential is enough for speci- 
fying the velocity field (See Eq. (13), Eq. (29) and the results 
shown in Sec. 3). In general, given the density distribution 
computed by the methods introduced in Section 2.2 together 
with the conclusion in Section 2.1, we are allowed to acquire 
the exact velocity field by Eq. (A. 17) in Casertano (1983). 
Fully solving the Poisson equation becomes not necessary. 
Fourth, initializing a disk over distributed memories allow 
us to deal with a larger data set which cannot be fully con- 
tained in a single memory storage. 

We propose a simple but very effective way of initializ- 
ing a three-dimensional gas disk. This method can be easily 
incorporated into any existing code based on either a La- 
grangian or Eulerian approach. No data exchange between 
processors is needed. Vertical density profile is obtained self- 
consistently without solving the full Possion equation. We 
implement these ideas with the adaptive mesh refinement 
magnetohydrodynamics code RAMSES (Teyssier 2002) and 
apply our concepts to probe the onset of the disk instability. 
We modify the dispersion relation for the infinitesimally thin 
disk (Lin &: Shu 1964) to be able to treat thick disks (Gol- 
dreich & Lynden-Bell 1965; Kim & Ostriker 2002a, 2006; 
Shetty & Ostriker 2006; Lisker & Fuchs 2009). The thresh- 
old value Qth is then obtained semi-analytically. Previous 
studies on this subject are either focused on a small patch 
of a galaxy (2D/3D: Kim & Ostriker 2002a) or are globally 
two-dimensional but with the reduction of gravity included 
in the governing equations (Shetty & Ostriker 2006). In this 
paper, we revisit the subject as a test of our fully three- 
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dimensional isolated galaxy models. Models with or without 
stellar potential are investigated. 

Galactic disks are comprised of stars and gas. Both com- 
ponents are coupled to each other via the Poisson equation. 
Since the stellar disk dominates the mass budget within the 
luminous disk, its presence has great impact on the scale- 
height of the gas disk as described in NJ02. A balanced 
initial condition depends not only on the correct vertical 
structure but also on the correct rotation velocity. To spec- 
ify the rotation velocity needed, the mass enclosed within a 
certain radius must be under control. Although it is com- 
mon practice to specify the functional form of the volume 
densities of 3D disks, we show that because of the flaring 
disk this typically results in a surface density profile that 
contains a central 'hole' (Agertz et al. 2009). This problem 
can be trivially avoided by specifying the desired surface 
density profile instead. We show that the corresponding vol- 
ume density can easily be obtained using a simple iterative 
scheme. The surface density of the total gas (HI + H2) from 
observation (Leroy et al. 2008) typically follows an expo- 
tential profile in disk galaxies. This profile gives an analytic 
description of the total mass enclosed within a radius as well 
as a reasonable approximation for velocity field as shown by 
Eq. (29) below (Binney & Tremaine 2008). 

Describing the stellar disk with a fixed background po- 
tential is at best an approximation to reality. The interaction 
between live stellar disk and gas can potentially destabilize 
the system (Rafikov 2001; Li, Mac Low & Klessen 2005a, 
2005b, 2006; Kim & Ostriker 2007). After all, the gas is cold 
compared to the stellar disk and has highly non-linear re- 
sponse to the asymmetric stellar potential. The gravitational 
interplay between the collisionless stars and dissipative gas 
is important for a number of key questions in galactic dy- 
namics. For example, what is the physical origin of grand de- 
sign spirals? Or what initiates and regulates the formation 
of stars? Having access to well-controlled initial and envi- 
ronmental conditions is a prerequisite to discovering their 
causes. 

This paper is organized as follows. The ideas of initial- 
izing a gas disk are outlined in Section 2. Details of the sim- 
ulation parameters and test runs are described in Section 3. 
Axisymmetric instability of the disk is revisited in Section 
4. The self-induced spirals due to swing amplification will 
be discussed in Section 5. A brief summary and the possible 
extension of this work is put in Section 6. 



2 FORMULATION OF EQUATIONS 

In this Section, we develop the required relations and equa- 
tions to immerse a 3D gas disk in a preexisting static poten- 
tial. Assuming that the gas disk and the preexisting poten- 
tial share the same symmetry axis, cylindrical coordinates, 
(R, 4>, z), are adopted to formulate the dynamics of the sys- 
tem. Axial-symmetry enables us to discard the terms de- 
scribing the variation in azimuthal direction, i.e., d/d(f> = 0. 
A gas disk which is in detailed balance should fulfill the 
following equations: 



R 



1 dp <9£ 
p g OR + OR 

p g dz dz 



(1) 
(2) 



with p g , p, Vrot and $ being the volume density of the gas, 
the gas pressure, the azimuthal rotation velocity ("rotation 
velocity" in short) and the joint total potential. Equation 
(1) describes that the gravitational pull in radial direction 
is counterbalanced by the centrifugal force and the pressure 
gradient. Equation (2) states that hydrostatic equilibrium 
along the symmetry axis, the ^-direction, is achieved by the 
balance between vertical pull of the gravity and the pressure 
gradient in z. 

To make the system self-consistent, the Poisson equa- 
tion must be involved: 



V 2 $ = 47rG(p g + pom + p s ), 



(3) 



with G, pdm and p s being the gravitational constant, and the 
volume density of dark matter and stars. The total potential 
is comprised of the contributions from the dark matter halo, 
the stellar disk and the self-gravity of the gas, i.e., $ = 
$dm + $ s + $g. In addition, the ideal gas law provides the 
link between the gas density, the gas temperature and the 
gas pressure: 



p = p g ( 7 -l)e(T), 



(4) 



where 7 represents the ratio of the heat capacities (adiabatic 
index) , e the specific internal energy and T the gas temper- 
ature. In the case of an ideal gas, the specific internal energy 
depends only on temperature, and is given by 



1 p,m p ' 



(5) 



with ItB being the Boltzmann's constant, p, the atomic 
weight and m p the mass of a proton. However, to close the 
set of equations, we should either invoke the energy equation 
or an equation of state (EoS), which will be used to evolve 
the system. 

A disk which is in hydrodynamic equilibrium should 
stay in its original state if we evolve the disk with the same 
equation of state which is used to set-up the disk. The 
numerical results throughout this paper are based on the 
isothermal equation of state, i.e., 



P = c s 



(6) 



with c a being the sound speed, a temporal and spatial con- 
stant. Equations (1) to (6) then form the basis of our discus- 
sion. In this paper, all the disks are in detailed equilibrium 
with the isothermal EoS. If those disks are adopted to evolve 
with a cooling function or a polytropic EoS, we can make 
sure any change in temperature or dynamics is purely caused 
by a cooling or a heating source. 

For a polytropic gas, p = Kp^. with F and K being 
constant, integrating Eq. (2) gives: 



p s (R,z) = p g (R,z = 0) 



r- 1 



d{R,z = Q) 



,(7) 



where <fr z (R,z) = &(R, z) — <b(R,z = 0) defines the verti- 
cal potential difference. We have used the fact that c 2 = 
dp/dp g = ATpg -1 when approaching Eq. (7). Note that 
given r^l, the internal energy has the following relation: 



e(T) 



Kpl 



7-1 



(8) 



Combining Eq. (5) and Eq. (8) gives the temperature field 
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as a function of position if the gas disk is initialized with a 
non-isothermal EoS. As a special case, when T — > 1, Eq. (7) 
then converges to a form for the isothermal gas: 

p s (R, z) = p s (R, z = 0) exp (- Z) ^ . (9) 

As we can see from Eq. (7) and Eq. (9), the vertical structure 
of gas disk depends on the gas properties in the midplane 
and the vertical potential difference. 

To fully characterize a gas disk which is in detailed bal- 
ance, we need to specify the velocity, the density and the 
temperature at every location in the beginning of the sim- 
ulation. In the following sub-sections we study the general 
properties of the velocity and density field, which allows us 
to devise a simple, but effective method to initialize a 3D 
gas disk in hydrostatic equilibrium. 



to be a very thin structure embedded in a static potential 
contributed by the background spherical dark matter and 
the stellar disk. Because the gas disk is observationally thin 
we neglect the radial variation compared to the vertical one 
(i.e., \(d/dR(Rd$ e /dR))/R\ < \d 2 $ g /dz 2 \).In Appendix E 
we show that this is a valid assumption for realistic gas disks. 
For an axisymmetric thin disk, the Poisson equation then 
reduces to (Binney & Tremaine 2008): 

^ = ^Gp g . (14) 

with $ g being the potentials contributed by the gas. In the 
following, we focus only on disks with initially constant tem- 
perature, i.e., the rotation velocity required for equilibrium 
has no dependence on the height above or below the mid- 
plane. 



2.1 Azimuthal Rotation Velocity 

In this sub-section, we treat the azimuthal rotation velocity 
as generally as possible. To make the notation concise, we 
drop the subscript of gas density, p g , and restore the sub- 
script after this sub-section. Without further assumption, 
integrating Eq. (2) from to 2 gives: 



f 

Jo 



ldp, 

~7T dz ' 

p dz 



By integrating Eq. (10) in parts, we have: 
p(R,z) _ p(R) 



p(R,z) p(R) 



r p_dp 

Jo P 29z 



dz-$ z (R,z). 



(10) 



(11) 



Inserting Eq.(ll) into (1), we get (see Appendix A for fur- 
ther details): 



vL{R,z) 



R 



1 dp 
pdR 



ldp_ 
pdR 

+ 

=0 

dp 



d<f>(R,z) 
OR 



dR 



rdp\ d_ (p_ 

\dz) dR \ p 2 



4 )Uz. (12) 



fdp\ d_ 

VdR) dz 



Equation (12) shows that the rotation velocity is indepen- 
dent of height above or below the midplane so long as the 
integral vanishes. It is evident that gas with a barotropic 
equation of state, i.e., p(p g ), fulfills this requirement. In ad- 
dition, for an initially constant temperature (T is everywhere 
the same in the beginning) disk, the initial pressure is a 
function of volume density only and therefore the integral 
becomes zero. In this case, equation (12) can be simplified 
further: 



V 2 ot (R,z) = R 



<9£ 
dR 



+ ( 7 -l)e 



d hip 



d\nR 



(13) 



Equation (13) states that the process of specifying the initial 
velocity in 3D comes down to knowing the rotation velocity 
in the equatorial midplane. 



2.2 Density Distribution 

From now on, to avoid confusion, we restore the subscript for 
the gas density. To proceed further, we consider the gas layer 



2.2.1 Density Method 

Differentiation Eq. (2) with respect to z and inserting 
Eq. (14) leads to a second-order non-linear differential equa- 
tion: 

d 2 p 1 d Pg dp ,. „ d 2 $ s d 2 $ r 



dz 2 p g dz dz 8V 8 



+ 



dz 2 



dz 2 



') = 0, (15) 



with <3> s and $dm being the potentials contributed by the 
stellar disk and the dark matter, respectively. So far, Eq. 
(15) is still general with respect to any kind of equation of 
state. However, a single equation with two unknowns p and 
p g is not solvable. To continue with Eq. (15), in this sub- 
section, we assume that the gas is barotropic, i.e., p(p g ). 
Given density distributions of stars, the dark matter and 
the boundary conditions in the midplane: 



Ps (R,0) =po(R) and ^ 
az 



0. 



(16) 



equation (15) can be solved by numerical integration, e.g., 
using the Runge-Kutta method. For a single- component, 
self-gravitating, locally isothermal disk (c B (R) can be a func- 
tion of radius), Eq. (15) has an exact solution: 



p s (R,z) = po(R)sech 2 (z/h) 



(17) 



with po{R) being the gas volume density in the midplane, 
h = cf /2ivG P o the scale-height and c s the local isothermal 
sound speed. According to Eq. (17) and since the midplane 
volume density, po(R), generally decreases with radius, to 
keep the scale-height a constant, the sound speed and there- 
fore the temperature must be a function of radius. 

Equation (15) is the simplified version of the formula 
derived by NJ02 (see also Kim et al. 2002a), where they in- 
vestigated the vertical structure in a gravitationally coupled, 
multi-component galactic disk. It is important to notice 
that all calculations can be done locally without the need 
of exchanging information between processors and therefore 
greatly reduces the complexity of coding. 

In principle Eq. (15) allows one to compute the den- 
sity of the gas such that the disk initially is in hydrostatic 
equilibrium. The actual implementation using Eq. (15) does 
not guarantee the positivity of the density. In particular, 
at large radii Pg (R,z) is typically close to zero, and small 
errors due to the numerical integration often yield negative 
densities. This problem is especially relevant when using the 
adaptive-mesh refinement techniques. 
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Initializing a gas disk with AMR usually starts with the 
coarsest grid. A natural selection of the integration step is 
the cell size. Then a problem immediately rises when solving 
Eq. (15) to specify the volume density. Supposing that the 
cell size is much larger than the scale-height of the gas disk, 
the errors introduced by the coarse integration may lead to 
negative densities on the outskirt of the computation do- 
main. One might think the integration can be done by using 
either adaptive integration intervals or simply a fixed inte- 
gration interval which is much smaller than the cell size. 
However, the improvements only guarantee the convergence 
of the solution not the positivity. Nevertheless, because of 
the generality of Eq. (15), density method is still valuable 
for semi-analytic research. 

2.2.2 Potential Method 

In this sub-section, we develop another route for specifying 
the density distribution. We stress that the following deriva- 
tion is only applicable to initially isothermal disks. With this 
constraint, integrating Eq. (2) gives: 



p g (R,z) = p (-R)exp 



( 7 -l)e 



(18) 



Combining Eq. (14) and Eq. (18), a second-order non- linear 
equation for the vertical potential difference of gas is ob- 
tained: 



d 2 $ E 



dz 2 



4irGpo(R) exp 



"(7-l)e 



(19) 



Given the analytic forms of $dm and <J> S the only unknown is 
the potential difference of gas, $ gi z = $ g (i?, z) — <fc g (R, z = 
0). Similar to the density method, given the boundary condi- 
tions po(R), <b(R, z = 0) and d$(_R, z = 0)/dz = 0, numeri- 
cal integration can be applied to solve Eq. (19). By inserting 
<E> Z obtained by integrating Eq. (19) into Eq. (18), the den- 
sity distribution is acquired. Notice that what really matters 
to us is the potential difference, not the absolute value. This 
means the value of $ g (R, z = 0) can be an arbitrary con- 
stant, although we do know the values of <&om(R, z = 0) 
and ®s(R,z = 0). 

The merit of this formulation is evident, the occurrence 
of negative density is avoided by Eq. (18). Tiny errors in the 
potential difference will not do any harm to the positivity 
of the gas density. Numerical experiments show that in nor- 
mal cases in which both the density method and potential 
method work, the solutions are consistent. 

At a given radius, R, solving Eq. (19) only provides us 
with information about the potential difference, &z(R, z). 
However, a useful byproduct of the potential method is that 
it is possible to acquire a good approximation of the total 
potential by: 



* S (R, z) = $ g (J2, 2 = 0) + * StZ {R, z), 



(20) 



as long as we know the potential in the midplane, & g (R, z = 
0). Equation (20) is an approximation because the use of 
Eq. (19) is based on the reduced Poisson equation Eq. (14) 
in which the variation in radial direction is ignored. The 
gradient of the potential <E> g (R, z — 0) determines the ve- 
locity field required while the vertical potential difference 
& gtZ (R, z) gives the vertical structure of the disk. In princi- 
ple, the radial force, which is associated with <fr g (R, z = 0), 



in the equatorial plane for an axially symmetric density dis- 
tribution can be evaluated precisely by the equation (A. 17) 
in Casertano (1983). This allows us to obtain the total po- 
tential without fully solving the Poisson equation. In prac- 
tice, if the initialization is performed with multi-node clus- 
ters, each node only keeps part of the information about the 
density distribution, data exchange with AMR itself is tech- 
nically challenging. In Section 3, for an exponential disk, the 
numerical results will show that the use of Eq. (29) is a good 
approximation for most of our interests. The corresponding 
$ g (R, z = 0) associated with Eq. (29) can be found in the 
book by Binney & Tremaine (2008), Eq. (1.164a). 

Equation (20) is useful, because involving the total po- 
tential into the formulation is an important step for self- 
consistently building up the combined disks comprised of 
a live stellar disk and a gas disk. Extension to the work 
of Shu (1969), Kuijkcn & Dubinski (1995, hereafter, KD95) 
develop a self-consistent disk-bulge-halo model for galaxies. 
The distribution function built by Eq. (6) in KD95 involves 
the potential differences <J> Z and <fr(R, 0) — &(R C , 0), with R c 
the radius of the guiding center. The potential method pre- 
sented here can be naturally incorporated into the frame- 
work of KD95. Therefore, in this paper, all the disks are 
initialized by the potential method. 



2.2.3 Exponential Disk 

Some studies have assumed that the midplane density of a 
3D gas disk has an exponential form (Tasker 2006, Agertz 
2009). However, as we now demonstrate, in general this re- 
sults in a surface density distribution that peaks at a specific 
non- zero radius, giving rise to a ring-like feature. We assume 
a gas disk with the popular sech 2 vertical profile: 



p s (R, z) = p c exp (—R/Rd) sech 2 



2 / Z 



h(R) 



(21) 



with p c being the central volume density, Rd the disk scale- 
length and h(R) the scale-height as a function of radius. The 
surface density then reads: 



T,(R)= p g (R,z)dz = 2p c exp(-R/R d )h(R). (22) 

J — oo 

Based on Eq. (22), we measure the scale-height of a disk 
at certain radius by h(R) = E(_R)/(2po(7?)). The extrema of 
the surface density can be evaluated by taking the derivative 
to Eq. (22): 

dg(g) „ , D/D , f dh(R) h(R)\ 

— = 2p c exp(-iVi? d ) ^_-_j=0. (23) 

We Suppose that the disk is linearly flaring, i.e., h(R) — h + 
R/Rh, with ho being the minimum scale- height of the disk 
and Rh a factor controlling the degree of flaring. The peak 
of the surface density then locates at i? poa k = Rd — hoRh- 
Whenever the i? p0 ak is positive, we get a ring in surface 
density. However, a ring in the surface gas density is not 
commonly seen in a real disk galaxy. An exponential pro- 
file in the total gas is prevalent in disk galaxies (Leroy et 
al. 2008). 

In order to avoid this feature, it is advantageous to spec- 
ify the actual surface density of the disk, rather than its 
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Table 1. Models' Parameters. 



Run 


T(K) 


M S (M Q ) 


Figure 


GasO 


4 x 10 4 




(1),(2),(3) 


Gasl 


2 x 10 4 




(4),(5) 


Gas2 


1 x 10 4 




(4),(5),(8) 


Gas3 


9 x 10 3 




(4),(5) 


Gas4 


8 x 10 3 




(4),(5) 


GasStarl 


7 x 10 3 


4 x 10 10 


(6),(7),(8) 


GasStar2 


6 x 10 3 


4 x 10 10 


(6),(7) 


GasStar3 


5 x 10 3 


4 x 10 10 


(6),(7) 


GasStar4 


4 x 10 3 


4 x 10 10 


(6),(7) 



"All disks have a gas mass of 10 10 A4q. 



midplane density. In the case of the exponential profile, the 
surface density reads: 



E(fi) = S exp(-i?/7? d ) 



./ — : 



p s (R, z)dz, 



(24) 



with E being the surface density in the galactic centre. 
Combining Eq. (24) and Eq. (18), the volume density in the 
midplane can be expressed as: 

E exp(- J R/i? d ) 



po(R) 



£>xp(-$ z /[( 7 -l)e])da 



(25) 



It shows that the correct volume density in the midplane 
for the desired surface density profile can be obtained iter- 
atively. Given a initial guess for po(R), & z is evaluated via 
Eq. (19) and also the integral appears in Eq. (25). One needs 
to iterate between Eq. (19) and Eq. (25). However, depend- 
ing on the quality of the initial guess, convergence can be 
reached very fast. For instance, with the initial guess being 
po(R) = Eo exp(— R/R d ), a six-time iteration already gives 
us a reasonable exponential disk. 

We pursue the exponential disk for several reasons. One 
is simply because it is commonly seen in disk galaxies. An- 
other is that we have a better control of the total mass. As 
we can see, if we specify the midplane volume density instead 
of the surface density, we do not exactly know the total mass 
until we finish the integration. Without knowing the total 
mass in advance, evaluating the circular velocity contributed 
by the self-gravity will not be a trivial task. Nevertheless, in 
principle, any profile of the surface density can be achieved 
simply by the process introduced in this sub-section. 



3 IMPLEMENTATION AND TESTS 
3.1 Simulation Parameters 

In this Section, we test the ideas outlined in the previous Sec- 
tion. We implement the method in the AMR-code RAMSES 
(Teyssier 2002). RAMSES uses grid-based Riemann-solvers 
for the magneto-hydrodynamics (MHD) and particle-mesh 
(PM) technique for the collisionless physics. It has a fully 
parallelized Poisson solver with periodic boundary condi- 
tions, which we use for this paper. Gas disks which are ini- 
tialized isothermally with an exponential surface density of 
a scale-length of 3.5 kpc and a total mass of 1O 1O M0 are 
embedded in a static potential. An isothermal equation of 
state is used to evolve the disks throughout this paper. 

The tests are mainly divided into two groups, one group 



is evolved with a static stellar potential (models with the 
prefix GasStar), the other without (models with the prefix 
Gas). Gasl to Gas4 are M33-like gas-rich galaxies, while 
GasStarl to GasStar4 are more similar to the Milky- Way. 
The main parameters of the models are listed in Table 1. 
The size of the computational domain is 250 kpc on a side. 
Up to 12 levels of refinement are used for those runs without 
stellar potential, and 13 levels for the other group, i.e., the 
corresponding highest spatial resolutions are about 60 pc 
and 30 pc, respectively. 

The volume density of the halo is described by the NFW 
profile (Navarro, Frenk & White 1997): 

I \ Mwo CX tna\ 
PDM(r) = -r^ 77— T77 , (26) 

47r/(c)r 20 o r 2 (l + a;) 2 ' 

with the Virial mass M 2 oo = 10 12 M Q , x = rc/r 2 oo, concen- 
tration parameter c = 12, distance r — \/R 2 + z 2 , Virial 
radius rzoo = 213 kpc and /(c) = ln(l + c) — c/(l + c). The 
Virial radius (r2oo) is a radius within which the averaged 
matter density is 200 times the critical density. 

The density distribution of the stellar disk reads 
(Miyamoto & Nagai 1975, Binney & Tremaine 2008): 



Ps(r) = 



b 2 M a \ aR 2 + (a + 3Vz 2 + b 2 ){a + Vz 2 + b 2 )' 



47T 



[R 2 + (a + Vz 2 + b 2 ) 2 ] 5/2 (z 2 + b 2 f/ 2 



,(27) 



with M s = 4 x 10 10 M Q being the mass of the stellar disk, 
a = 3.5 kpc and b = 0.2 kpc the shape parameters. 

In light of the result drawn from Section 2.1, for an ini- 
tially constant temperature setup, we only need to know the 
circular velocity in the midplane for initializing the velocity 
field. The rotation velocity, V TO t, is decomposed into four 
components: 



V 2 

^rot 



vg M + V 2 + V 2 + V 2 



(28) 



where Vdm, V b , V g are the circular velocities corresponding 
to the dark matter halo, the stellar disk and the gas disk, and 
V p denotes the contribution due to the pressure gradient. 

In this paper, we have the analytic form for Vdm and V s . 
For the contribution from the gas disk and pressure gradient, 
we take the approximation for an infinitesimally thin disk 
with exponential surface density as described in Eq. (24). 
We set: 



V 2 (R) = ^G^ R d y 2 [I {y)K {y) - h^K^y)] 



V P \R) = ( 7 -l)e 



9 hip 



dlnR 



(29) 
(30) 



with y = R/(2Rd), Io, Ko, I\ and K\ being the modified 
Bessel functions of the first and second kinds of zeroth/first- 
order, respectively. Equation (30) derives from the second 
term of Eq. (13). However, contribution from pressure gra- 
dient in the midplane can only be evaluated after the gas 
disk is set up. Note that for an exponential disk, surface 
and volume densities decrease with radius and hence V 2 is 
negative. 



3.2 A Stable Disk 

To demonstrate that the disk built by the potential method 
described in Section 2 is in detailed equilibrium, we start 
with a stable equilibrium disk in model GasO. In this test, the 
stellar disk is deliberately removed. Without the dynamical 
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5 10 15 5 10 15 

R [Kpc] R [Kpc] 

(a) (b) 



Figure 1. Model GasO: (a) The total rotation velocity (solid) and contributions from dark matter halo (dashed), gas (dotted) and gas 
pressure (dash-dotted). Note that we plot the absolute value of the pressure gradient to have positive values for the direct comparison. It 
should be in opposite sense to the gravity. In this model contributions from the gas self-gravity and the pressure gradient is not negligible, 
(b) The Q value of model GasO as a function of radius as defined by Eq. (32). The Q is well above the threshold value Q t h = 1, thus the 
disk is expected to be stable. No structure should develop with time. 



1.601 Gyr 




5 10 15 5 10 15 

R [Kpc] R [Kpc] 

(a) (b) 

Figure 2. Model GasO: The evolution of 1.6 Gyrs of a stable disk, (a) The evolution of the surface density (upper panel) and the 
rotation curve (lower panel) at t = Gyr (solid) and t = 1.6 Gyr (diamond). Overall, the surface density and rotation curve are kept 
very well over 4 orbital periods, (b) The evolution of the scale-height at t = Gyr (solid) and t = 1.6 Gyr (dotted). The small change in 
scale-height indicates that the required circular velocity is overestimated probably due to the approximation of Eq. (14) and Eq. (30). In 
all, the disk still stays well in the initial condition. The step-wise character of the scale-height reflects our discretization and the change 
of spatial resolution due to the AMR. 



support from the stellar disk, the self-gravity of the gas plays 
the dominant role in determining the vertical structure of 
the disk and provides a not negligible contribution to the 
rotation velocity. 

Figure la decomposes the rotation curve into the dif- 
ferent contributions by the halo, the gas and the pressure 
gradient. Note that the forces of the self-gravity and the 
pressure gradient are in opposite sense, the self-gravity pulls 
matter inwards while the pressure gradient pushes outwards. 
In this figure, V p is shown in its absolute value. If we ignore 
the pressure gradient, the disk would rotate too fast and 
gradually drift outward. Figure lb shows the conventional 
Toomre's Q defined by: 

« = iSt (31) 

with ft being the epicyclic frequency. It shows that the Q is 
well above Qth = 1, the threshold value for stability, at all 



radii. The disk is hot enough to keep the disk stable and no 
structure should develop with time. 

We let the disk evolve for 1.6 Gyrs (four orbits for the 
gas at 10 kpc) and check how well the disk properties are 
kept. Figure 2a presents the evolution of the surface den- 
sity and the rotation curve. The solid lines represent the 
initial states and the diamond symbols are the status after 
an evolution of 1.6 Gyrs. The surface density is obtained 
by projecting along the symmetry axis and the rotation 
curve is evaluated by the mass-weighted circular velocity, 
v T ot(R) = J Pg(R, z)v vo t(R, z)dzfSg(R). Although a small 
amount of mass accretes onto the very central part of the 
disk, overall the surface density and the rotation curve are 
kept very well. Mass accretion into the center seems unavoid- 
able for a Cartesian-grid code mainly because too small a 
number of cubic cells is used to mimic the circular motion 
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in the centre. This accretion will be eventually halted by the 
pressure gradient built by the accumulating material. 

Figure 2b shows the evolution of the scale-height. The 
solid line represents the initial state and the dotted line 
the evolution after 1.6 Gyrs. Upon closer inspection we 
find that the disk undergoes a very small amount of mass- 
redistribution in the radial direction, which we believe to 
be a consequence of our two approximations when initializ- 
ing the disk. One is from the reduced Poisson equation, Eq. 
(14), and the other is from the use of Eq. (29). Equation 
(29) overestimates the circular velocity needed to support 
the disk. The thickness of the disk reduces the potential in 
the midplane by a few percent (see appendix B). Figure 3 
shows the snapshots of the face-on surface density at t = 
(Fig. 3a) and at t — 1.6 Gyr (Fig. 3b). No structure is de- 
veloping during the course of the simulation. 

To sum up, figures 1 to 3 indicate that without external 
perturbation the disk is quiet over secular time-scales. The 
shape of such a disk is naturally flaring, i.e., the scale-height 
increases with radius. The ideas described in Section 2 are 
able to treat the initial condition self-consistently. A well- 
balanced disk is especially useful to probe the onset of disk 
instability as described in the next Section. 



4 AXISYMMETRIC INSTABILITY 

The question of disk stability has been investigated for more 
than four decades since the pioneering works by Toomre 

(1964) for collisionless stars and Goldreich & Lynden-Bell 

(1965) for gas sheets. Understanding the origin and evolu- 
tion of disk structure is challenging. If the disk is stable 
like our model GasO, no structures can form. On the other 
hand, if the disk is highly unstable, the surface density will 
quickly fragment and develop a clumpy and chaotic-looking 
appearance. There will be no well-organized structures. The 
striking spiral appearance of many nearby disk galaxies in- 
dicates that those disks are marginally stable. 

For an infinitesimally thin disk, the instability thresh- 
old is at Qth = 1 (Toomre 1964). The first theoretical work 
to include the finite thickness of a self-gravitating gas disk is 



that by Goldreich & Lynden-Bell (1965). Some authors have 
investigated the stability of finite thickness gas disks in nu- 
merical simulations (both in 2D and 3D) using local patches 
within a shearing box (Kim & Ostriker 2006, 2002a; Gam- 
mie 2001). This technique, in 2D, has also been used by Kim 
& Ostriker (2007) to investigate the interaction between the 
gas disk and a live stellar disk. Shetty & Ostriker (2006) used 
global 2D simulations in which they incorporated the effect 
of finite disk thickness by diluting the gravitational force. 
For 3D global disk calculations, see Li, Mac Low & Klessen 
(2005a, 2005b, 2006), who investigate the relation between 
disk instability and star formation rate. These studies all 
agree that although the inclusion of the thickness does not 
have a qualitative impact on the disk instability, it does shift 
the threshold value of instability quantitatively. In addition, 
accounting for disk thickness may have a large impact on 
the evolution of a disk, such as the development of spurs or 
the wiggle instability (Kim & Ostriker 2002b, 2006). 

In this paper, armed with a well-balanced gas disk, we 
revisit the axisymmetric instability of disks in 3D global 
fashion. We first derive the reduction factor F which re- 
flects the reduction of the gravity due to the finite thickness 
of the disk. Then the corresponding instability threshold 
Qth(R,T) derived from a semi-analytic calculation is com- 
pared with the numerical results. In the final sub-section, 
we also explore the impact of the presence of a static stellar 
potential on the axisymmetric instability. 

4.1 The impact of thickness on disk stability 

The Fourier component of the perturbed gravitational po- 
tential, $fc, of an infinitesimally thin disk is given by: 



where k represents the wave number of the Fourier compo- 
nents and x — R — Rq being the radial deviation for an ax- 
isymmetric perturbation. Supposing that a 3D disk is piled 
up by a stack of infinitesimally thin gas layers, we approx- 
imate the effect of the disk thickness by superimposing the 
contribution from every razor-thin layer: 
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27rGS fc e'- 
1*1 



./ — C 



„-k\z-h\ sech 2 (h/ h z ) Alt 
e 2h z dh ' 



(33) 



with h z being the scale- height of the disk. In Eq. (33), we 
model the vertical structure of the gas disk by a sech 2 func- 
tion. This is valid especially for the inner part of disks 
where the vertical structure is mainly determined by the 
self-gravity of the gas. See also the Fig. Dl in Appendix D. 
Equation (33) leads to the Fourier potential in the midplanc: 



$fc(z = 0) 



27rGS fc e' 
\k\ 



-F(k,h z ), 



(34) 



with F(k,h z ) being the reduction factor described by (see 
Appendix C): 



with H being the harmonic number defined by: 



F(k,h z ) = 1- -kh 



H(a) 



1 



V 



-Ay. 



(35) 



(36) 



The Lin-Shu (1964) dispersion relation for the axisymmetric 
perturbation is then modified to: 



lo 2 = n 2 — 2nCEo\k\F(k, h z ) + c 2 k 2 . 



(37) 



The dispersion relation states that on small scales (k — > oo) 
the disk is stabilized by gas pressure, i.e., the term c 2 k 2 . 
Large scales (k — > 0) are regulated by global shear, i.e., the 
k 2 term. The instability however happens at intermediate 
wavelengths, much smaller than the disk size but still larger 
than the thickness of the disk. In this region, neither global 
shear nor gas pressure can resist the gravitational collapse. 
The reduction factor, < F < 1, softens the effect of self- 
gravity and makes the disk more stable. 

Given a certain radius R and temperature T, we ob- 
tain the threshold value Q t h(T, R) by probing the maximum 
value along the neutral curve defined by setting uj 2 = in 
Eq. (37) and calculating the epicyclic frequency, k, from 
the rotation curve. Similar to the conventional Toomre cri- 
terion for the stability of an infinitesimally thin disk, Q t h 
is a threshold curve for thick disks. Above Qth the disk is 
stable and otherwise unstable. Since the Qth is a function 
of both temperature and radius, it is convenient to define 
the critical value Q cr it, which is the value of Q th for which 
Qth(T,R)/Q(R) — 1, and the corresponding critical tem- 
perature Tcrit- 

The solid lines shown in Fig. 4 represent the threshold 
value Qth as a function of radius. Each plot corresponds to a 
disk of different temperature. The dash-dotted lines are the 
actual Q values defined by Eq. (31) of the different models. 
From these figures, the most unstable radius is around R = 2 
kpc. The corresponding surface densities after an evolution 
of 750 Myrs are shown in Fig. 5. The gas at the most unsta- 
ble region has revolved for more than four orbital periods 
around the disk center. 

These figures shows that the prediction of Q cr it and the 
numerical results match quite well. The Q value of Gasl 
is well above the solid line and shows a featureless surface 
density. As shown in Gas2 and Gas3, with the decrease in 
temperature, the Q th curves shift up and the disks' Q curves 
come down. As a consequence, the disk starts to develop 
multi-armed structure, which is very likely caused by swing 



amplification, as discussed in Section 5. And finally in Gas4, 
the curves Q and Qth intersect. The disk fragments and 
starts to behave chaotically. A more detail calculation shows 
that the two curves just touch each other at a temperature 
Tcrit = 8.5 x 10 3 K with the maximum threshold Q cr it = 
0.693, which is close to Q cr it = 0.676 of Goldreich & Lynden- 
Bell's (1965) analysis but away from the numerical result, 
Qcrit = 0.647, of Kim et al. (2002a). However, the actual 
value of Qcrit is model dependent. Different models of the 
dark matter, the stellar disk and even the EoS will all affect 
the resulting value of Q cr it- 



4.2 The inclusion of stellar potentials 

The inclusion of a static stellar disk alters two important 
factors which influence the stability of the disk. One is the 
rotation curve and the other is the thickness of the gas 
disk. By changing the rotation curve, the epicyclic frequency, 
k, changes accordingly. Supposing a flat rotation curve de- 
scribed by Q = Vo/R, the epicyclic frequency k then reads: 



2V 2 



2YI 



(38) 



with Vo being the rotation velocity. The presence of a stel- 
lar disk tends to stabilize the gas disk via increasing Vo- 
However, by increasing the gravitational pull in the vertical 
direction, the gas disk becomes thinner and therefore more 
susceptible to gravitational collapse. In Section 4.1, we have 
already seen that the scale-height, which is governed by the 
temperature of the disk, is a very sensitive factor for the 
disk stability. GasStarl to GasStar4 are designed to explore 
the competition between the two opponents. 

From Fig. 6, we first notice that, compared to Fig. 4, the 
threshold value, Q cr it, is boosted from 0.693 to 0.75 due to 
the decrease in scale-height. This makes the disk more prone 
to gravitational instability. On the other hand, the change 
of the rotation curve drastically shifts the dash-dotted curve 
upwards. Instability only sets in once the temperature of the 
gas disk drops below T cr it ~ 6000 K. Overall, the presence 
of the static stellar disk tends to stabilize the disks. 

Figure 7 shows the surface density after an evolution of 
250 Myr. During this period, the gas in the most unstable 
region has finished 2.5 orbits. All the gas disks are devel- 
oping multi-armed spiral structures within the region where 
the disk is the most vulnerable to instability according to 
Fig. 6. At this moment, the most unstable disk, GasStar4, is 
experiencing fragmentation. High density filaments are evi- 
dent from the image. While GasStar2 is still in its early stage 
of instability, GasStar3 is just about to enter the fragmenta- 
tion phase. GasStarl, on the other hand, does not fragments 
at all during the course of simulation. 

The trend is clear. The cooler the disk, the faster it frag- 
ments. The spiral structure seen in these images are due to 
swing amplification (Toomre 1981; Goldreich & Lynden-Bell 
1965), a mechanism that is capable of amplifying the per- 
turbation by swinging the leading waves to trailing. Swing 
amplification is effective as the disk Q (dash-dotted line) 
is approaching the threshold Qth (the solid line). The spi- 
rals are sheared, become tighter and tighter and enhanced. 
Once the density reaches the supercritical point, instability 
sets in. 
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Figure 4. Plots (a) to (d) correspond to models from Gasl to Gas4, respectively. In each plot, curves of the disk Q (dash-dotted) and 
the threshold value Q th (solid) are put together to probe the onset of axisymmetric instability. Q t h(-R) is a obtained by probing the 
maximum value along the neutral curve for a given radius. Information of the disk thickness has been encapsulated in the reduction 
factor defined by Eq. (36). When the two curves meet, we expect the disk fragments very fast. This figure shows that the most unstable 
region is about the radius R = 2 kpc. The fact that the Q t h curves are well below unity shows the impact of the disk thickness on the 
disk stability. 



5 SPONTANEOUSLY INDUCED SPIRAL 
STRUCTURE 

An interesting feature which is hard to ignore in Fig. 5 
and Fig. 7 is that the marginally stable disks are sponta- 
neously developing multi-arm spiral structures. We have al- 
ready seen in Section 3 and 4 that the effect of the disk 
thickness is to shift the range of the marginally stable re- 
gion downwards and therefore to stabilize the disk. As we 
systematically lower the temperature to probe the onset of 
instability, runs with as well as without stellar potential are 
experiencing swing amplification. 

Hohl (1971) found that disks which are marginally sta- 
ble to axisymmetric perturbation are prone to develop a 
large-scale bar structure. This finding initiated both numer- 
ical (Zang & Hohl 1978; Sellwood 1981, 1985; Fuchs & von 
Linden 1998; Sellwood & Moore 1999) and theoretical stud- 
ies (Kalnajs 1978; Sawamura 1988; Vauterin & Dejonghe 
1996; Pichon & Cannon 1997; Evans & Read 1998; Fuchs 
2001) of marginally stable disks. Goldreich & Lynden-Bell 
(1965) and Toomre (1981) pointed out that self-gravitating, 
differentially rotating disks are able to amplify spiral waves 
by shearing a leading wave into a trailing one. Three key in- 



gredients, self-gravity, shearing and epicyclic motions work 
harmonically to make the phenomenon now coined with the 
name 'swing amplification' happen. 

Three necessary conditions need to be fulfilled in order 
to facilitate the swing amplification (Toomre 1981; Fuchs 
2001; Fuchs & von Linden 1998; Binney & Tremaine 2008). 
First, the disk must be marginally stable, i.e., for an in- 
finitesimally thin disk, 1 < Q < 2, as defined by Eq. (31). 
Second, the parameter X = k CI i t R/m = fc C rit/fc y (Toomre 
1981; Binney & Tremaine 2008), with m being the num- 
ber of arms and fc C rit = K 2 /(27rGE g ) the critical wave num- 
ber, has to be of order unity, i.e., somewhere between 1 and 
3 (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966; 
Toomre 1981). Third, there must be a mechanism that is 
able to induce leading arms in the system either explicitly 
by hand (Toomre 1981) or implicitly by random fluctuation 
induced by numerical noise (Toomre 1990; Sellwood & Carl- 
berg 1984; Fuchs 2001). We notice that most of these works 
mentioned above are for live stellar disks not directly for the 
gas disk. But since the amplification principles are the same, 
the results are still applicable to pure gas disks. 

As shown in Fig. 8a and 8b, GasStarl gets more arms 
than Gas2 does. To be more quantitatively, Fig. 8c and 
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Figure 5. Images (a) to (d) correspond to models Gasl to Gas4, respectively. They show the face-on surface density at t = 750 Myr. 
The size of the images are 20 kpc X 20 kpc. The gas at the most unstable radius has orbited around the center for more than four 
times, (a) Since the disk Q is well above the threshold value Q t h> the disk is featureless. In models Gas2(b) and Gas3(c) the disk Q is 
approaching Q t h around R = 2 kpc, both disks are developing self-induced spirals due to swing amplification, (d) The disk fragments 
very fast once Q and Q t ^ intersect. 



8d show the Fourier components as a function of radius. 
They are obtained by doing Fourier transform to (E g (Ji, <j>) — 
Eg(i?))/E g (i?), where E g (J?) is the averaged surface density 
of a given radius. Note that the dominating modes tends 
to be multiples of m = 4. This is a consequence of using 
a Cartesian grid, for which m — 4 is the natural mode. 
However, the dominating mode is determined by physics. 
The dominating mode of Gas2 is m = 8 while in GasStarl 
m = 12. 

As is apparent from Eq. (38), including a stellar disk 
causes an increase in fc C rit ■ Consequently, a larger value of m 
is required in order to satisfy 1 < X < 3. From the image 
shown in Fig. 8a and the relation, fc cr it oc k 2 , to keep X 
a constant, the number of spiral arms in GasStarl can be 
crudely estimated as m ~ 15. More precisely, the number of 
spiral arms, m, is predicted by (Toomre 1981; Athanassoula, 
Bosma & Papaioannou 1987; Fuchs 2001, 2008): 



with A n 



2-kR 



x being defined by: 

Acrit 



x (A/ny 



(40) 



(39) 



where A cr it = 27r/fc cr i t . The coefficient \ is a function of 
rotation curve (Fuchs 2001), as measured by Oort's constant 
A. 

We employ Eq. (39) to analytically estimate the number 
of arms and compare the predictions with the images shown 
in Fig. 8. For Gas2, spirals appear between 2 kpc and 5 
kpc. Within this radial range, the most unstable wavelength 
ranges from 2.0 to 3.6 kpc. The corresponding prediction 
for m ranges from 6 to 9, while the simulation reveals a 
spiral pattern with 8-fold symmetry. For GasStarl, spirals 
are prominent between 3 and 4 kpc, while the corresponding 
most unstable wavelength ranges from 1.4 to 2.0 kpc. The 
twelve arms developing in GasStarl should be compared to a 
predicted m ranging from 13 to 14. Hence, overall the trends 
in the simulations are in reasonable agreement with our pre- 
dictions. Note that the spatial resolution in both simulations 
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Figure 6. Plots (a) to (d) correspond to models GasStarl to GasStar4, respectively.: The Q (dash-dotted) and Q t ^ (solid) curves of 
the gas disks of different temperatures. The presence the stellar potential stabilizes the disks through changing the rotation curve and 
destabilizing the disk by increasing local gravitational force. The effect of disk thickness is included via the reduction factor Eq. (36). 
We need to lower the temperature down to T= 7 X 10 3 K in order to probe the onset of axisymmetric instability. Overall, the presence 
of the stellar potential stabilizes the disk. 



ranges from 60 pc to 120 pc, indicating that the most un- 
stable wavelengths are well-resolved. 

The observed small deviations can be explained as fol- 
lows. First, the formulation used to predict the number 
of arms is precise only for stellar disks. However, Toomre 
(1981) has shown the strikingly similar behavior of gaseous 
disks (Goldreich & Lynden-Bell 1965) and stellar disks (Ju- 
lian & Toomre 1966). Therefore, we have confidence that 
Eq. (39) is still applicable to gaseous disks. Second, the 
number of arms has to be an integer, a number of fraction 
given by Eq. (39) has no physical meaning. Third, the usage 
of a Cartesian grid introduces the multiples of the natural 
m — 4 mode, which manifests itself in the Fourier transform 
of the surface density. Fourth, swing amplification picks up 
the dominating mode. It takes some time to fully develop 
the dominating mode. All these factors combined determine 
the number of spiral seen in our simulations. It is impor- 
tant to realize that the most unstable radius according to 
the axisymmetric instability criterion might not be the most 
effective site for swing amplification, since the shear plays 
an important role in this process. 

Without any external pumping source, spiral waves pro- 
duced by swing amplification should be a transient phe- 



nomenon. Similar to material spirals, swing amplified spi- 
ral waves will experience azimuthal shearing which reduces 
their pitch angle until they become too tightly wound to be 
identified. As an example, in the Gas2 simulation, the spiral 
arm that appears around R — 2 kpc initially has a pitch 
angle of 90° and should be sheared to less than 1° within 
2.2 Gyr. On the contrary, we find that the spontaneously in- 
duced spirals seen in Gas2 can last for more than 3 Gyr and 
still keep the pitch angle relatively open. This result sug- 
gests at least one mechanism keeps replenishing noise into 
the disk, leaving the physics to pick up the dominating mode 
and sustain the waves. This noise can be caused by numerics 
or preexisting waves. 



6 SUMMARY 

In this paper we have developed a simple and effective 
method to compute the three-dimensional density and ve- 
locity structure of an isothermal gas disk in hydrodynamic 
equilibrium in the presence of an arbitrary external potential 
(i.e., dark matter halo and/or stellar disk). This is ideally 
suited to set-up the initial conditions of a three-dimensional 
gas disk in equilibrium in hydrodynamical simulations. We 
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Figure 7. Images (a) to (d) correspond to models GasStarl to GasStar4. They show the face-on surface density at t = 250 Myr. The 
size of the images are 20 kpc X 20 kpc. The gas at the most unstable radius has orbited around the center about two and half times. 
Spirals seen in model GasStar2(b) and GasStar3(c) are due to swing amplification. In (d) the disk fragments very fast mainly due to 
both the axisymmetric instability and swing amplification. 



first notice that as long as the gas is barotropic or has con- 
stant temperature at t = 0, the circular velocity needed to 
support the self-consistent disk is independent of the height 
above or below the midplane. This feature greatly simpli- 
fies the process of specifying the initial velocity field. All we 
need to know is the rotation velocity in the midplane. 

To specify the density distribution self-consistently, the 
hydrostatic equation coupled with the reduced Poisson equa- 
tion is adopted to develop the vertical structure of the gas. 
Two sets of second-order non-linear differential equations 
are found. One is directly associated with the gas density 
called the density method, the other associated with the gas 
potential called the potential method. In a simulation in- 
volving a huge dynamic range (using AMR techniques), the 
potential method is shown to be numerically more stable. 
A simple local iteration can be performed to gain a better 
control on the shape and the mass of disks. These ideas are 
simple enough to be incorporated into any existing code, 
and most importantly they are very effective. 

With gas disks that are in detailed balance, we are able 
to systematically investigate the axisymmetric stability of a 



fully three-dimensional disk for the first time. We probe the 
onset of instability both semi-analytically and numerically. 
Simulations without stellar disk show that the thickness of 
the gas disk, which is governed by the temperature of the 
disk, has a huge impact on the disk stability. The reduc- 
tion of the gravity decreases the threshold value by around 
30 percent in our models. As we gradually lower the gas 
temperature, the threshold Qth shifts up, the disk Q shifts 
down, and the system starts to develop multi-arm structure 
via swing amplification. The onset of the instability in simu- 
lations matches the theoretical prediction very well as shown 
in Fig. 4 and Fig. 5. The disk fragments as the two curves, 
Q and Qth, come very close to each other. 

The influence of the stellar disk is less obvious. Its pres- 
ence has a stabilizing effect on the gas disk through changing 
the rotation curve and a destabilizing one through the in- 
crease of the local gravitational force. The simulation results 
show that overall the presence of the stellar disk tends to sta- 
bilize the gas disk. But this conclusion comes with a caveat. 
The interaction between live stars and gas might be impor- 
tant. A live stellar disk itself can be unstable or marginally 



14 Wang et al. 




'°g, n <M c ../kpc^ 



7.5 



ll 6 - 5 




(a) Gas2 



lb) GasStarl 




4 6 
R [Kpc] 

(c) Gas2 



3.1 


25 


0.08 


20 


0.06 


E 15 


0.04 


10 


0.02 


5 








0.15 




0.05 



4 6 
R [Kpc] 

(d) GasSlarl 



Figure 8. The image size in (a) and (b) is 20 kpc X 20 kpc. (a) The surface density of Gas2 at t = 750 Myr. (b) The surface density of 
GasStarl at t = 500 Myr. In both cases, the inner parts of the gas disks, which have been evolved for about four orbital times, developing 
spiral structure. Contour plots (c) and (d) are the Fourier maps of (a) and (d), respectively. In (c) and (d), the horizontal axis represents 
radius, the vertical axis is the number of arms, m, obtained by Fourier analysis. The color represents the intensity of each Fourier mode, 
the redder the stronger. 



stable. Perturbations from the interstellar medium can trig- 
ger instabilities in the stellar disk. Since stars dominate the 
mass budget in Milky-Way type galaxies (more than 90 per- 
cent), and because gas is highly responsive and dissipative, 
the interplay between both components is one of the most 
interesting subjects in galactic dynamics. Tackling this prob- 
lem needs elaborate initial conditions for the live stellar disk 
or the combined disk. We stress that the potential method 
developed in this paper is compatible with the formulation 
in KD95. This makes the self-consistent combined disk a 
natural direction for future work. 

Marginally stable disks are susceptible to the process 
of swing amplification, a prevalent mechanism that triggers 
self-induced spirals. Simulations Gas2 and GasStarl show 
the spirals are prominent in the regions in which the gas 
can respond to swing amplification. Semi-analytic result re- 
lates the most vulnerable wavelength in azimuthal direction, 
Amax, to the number of arms. Numerically, The natural 
mode of a Cartesian grid together with the swing amplifica- 
tion determine the dominating mode of the spiral structure. 
Our numerical results with or without stellar disk shows the 
correct characteristics of the swing amplification. It happens 
in marginally stable disks and the number of arm fits reason- 
ably well to the analytic prediction. In the run of GasStar2, 



swing amplification eventually leads to disk fragmentation 
once the density becomes supercritical to the gravitational 
instability. However, in a sub-critical case like Gas2, the spi- 
ral structure can survive more than 3 Gyrs without frag- 
menting the disk, suggesting at least one mechansim is sus- 
taining the waves. The number of arms suggests a charac- 
teristic wavelength relating to the upper limit of the mass 
of giant molecular clouds (Escala 2008). 
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APPENDIX A: THE DERIVATION OF 
ROTATION VELOCITY 



Equation (11) can be re-written as 
p(R) 



p{R,z) = p(R,z 



p{R) 



- p(R,z] 

z=0 Jo 

p(R,z)[*{R,z)-*(R,z = 0)].. 



P dp 

p 2 dz 



dz 



(Al) 



where we have replaced $ z = &(R, z) — &(R, z = 0). Insert- 
ing Eq. (Al) in Eq. (1) involves a partial derivative to the 
integral, let us prepare it first: 
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With Eq. (A2), the first term of Eq. (1) then becomes: 
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Equation (11) says that the term in the big brace should 
vanish. And therefore, Eq. (1) reduces to 
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For the barotropic gas, i.e., p(p), the integrand of the integral 
vanishes: 
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dp 
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= 0(A5) 



For the cases of initially constant temperature, the specific 
internal energy, e, is a constant and therefore the pressure 
is a function of density only, the integrand vanishes. 



APPENDIX B: THE EFFECT OF THE DISK 
THICKNESS ON THE MIDPLANE POTENTIAL 

For an axisymmctrically and infinitesimally thin disk, the 
potential can be evaluated by the following relation (Binney 
& Tremaine, 2008): 



*(R,z)= [ 
Jo 



dkSo(k)J (kR)e 



-k\z\ 



(Bl) 



where Jo is the Bessel function of the first kind of order zero 
and So is the Hankel transform of —2nCSo defined by: 



So(k) 
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With Eq. (Bf) and Eq. (B2), we can superimpose the poten- 
tial contributed by each gas layer. For the sake of simplicity, 
we assume that the volume density has the double exponen- 
tial profile: 



p (R,z) = X e- R/R *- 



z/h z 



2h z 



(B3) 



with h z being the scale-height of the gaseous disk, Eq. (B2) 
then becomes (Gradshteyn & Ryzhik f 965, hereafter GR65, 
6.623-2): 
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with £ = 1/Rd- Az represents the infinitesimal thickness in- 
troduced to keep the dimension correct. The potential which 
takes into account the thickness of the disk then reads: 
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Evaluating the potential at the midplane, z = 0, yields: 
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Figure Dl. The force ratio F z ^yi/F ZtSZS at R=2, 5 and 8 kpc. 
It shows that the vertical structure of the inner disk is determined 
mostly by the self-gravity of gas. 

APPENDIX C: THE DERIVATION OF THE 
REDUCTION FACTOR 

To derive the reduction factor F defined by Eq. (35) we need 
to evaluate the integral of the form: 



-k\h\ c 



sech 2 (ah)dh = 2 / e" fe,l sech 2 (a/i)dft 



The last line can be reached by looking up the formulae 
3.541, 8.370, 8.361-7 listed in the integral table (GR65) and 
the definition Eq. (36). In the last line, we have employed 
the recursive relation (8.365-1 GR65): 



H(a) = H(a-l) + 



(C2) 



The asymptotic behavior of the harmonic number reads 
(8.367-2, 8.367-13 GR65): 



H{a) = lna + 7+ ia 1 



+ To7T Q + O(o? ),(C3) 



12 120 

with 7 = 0.5772156649 (8.367-1 GR65) being the Euler- 
Mascheroni constant. Note that Eq. (C3) is only reliable 
when a > 1. We employ the recursive relation (C2) to eval- 
uate H(a) for — 1 < a < 1. 



APPENDIX D: THE VERTICAL FORCE RATIO 

The vertical force ratio measures the impact of the halo force 
on the vertical structure. The simplified Poisson equation for 
isothermally self-gravitating gas disk reads: 
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Given the finite scale-height, the integral can be eval- 
uated numerically and compared with the result of the in- 
finitesimally thin disk. 



where h z being a measure of the scale- height. Parameter h z 
can be related to the volume density in the midplane, po(R) 
by: 
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Figure El. Contour map of e. The black lines represent the 
scale-height of the gas disk. 



The corresponding vertical force for the gas then becomes: 

F z>gas = —s- = -4nGh g po(R)tanh(z/h z ). (D3) 
oz 

For a NFW halo, the vertical force can be written down 
directly: 



F 2 



GM200 / c y x/(l + x) - ln(l + ; 
/(c) \r 2 oo 



,(D4) 



with x — c\/R 2 + z 2 /r2oo- Figure Dl then shows the force 
ratio F Zi om/ '-Fz, g as as a function of vertical height \z\ at dif- 
ferent radii. Comparing to the rotation curve shown in the 
left panel of Fig. 1, although the dynamics is still dictated 
by the potential of the dark halo, the vertical structure of 
the gaseous disk is mainly determined by the self-gravity 
of the gas component. However, it raises another issue, the 
presence of the stellar disk will dominate both the dynam- 
ics and the vertical structure of the gas and will affect the 
stability of the gas component via changing the thickness of 
the gaseous disk and the rotation curve. 



Here a is a constant that controls the scale-length of the 
disk and b(R), which we take to be a function of radius, 
modulates the scale-height of the disk. In the limit b — > 
this model reduces to the infinitesimal Kuzmin disk (e.g., 
Binney & Tremaine 2008). In an attempt to model the gas 
disk in our simulation 'GasO', we adopt a = 3.5 kpc. In 
order to mimic the flaring of the GasO disk (see Fig. 2b), we 
consider 



b{R) 



-1.58 x lO^T? 4 + 1.21 x 10 



~ 2 R 2 



0.20. 



(E3) 



Using the Poisson equation to solve (numerically) for 
the corresponding density distribution yields the radial- 
dependent scale-height shown as the solid black lines in 
Fig. El, and which is comparable to that of the GasO disk. 
The contours in Fig. El are defined by constant values of e. 
These show that our assumption that e 1 is well-justified 
in the inner part of the disk, out to ~ 3 scale-lengths, which 
encloses most of the disk mass. The assumption that e C 1 
deteriorates at larger radii and at higher altitude away from 
the midplane. This might be in part responsible for the very 
slight outward drifting of the disk seen in Fig. 2b. In cases 
that include a stellar potential and/or cooler gas, the gas 
disk is even thinner than the case considered here, result- 
ing in values for e that are even smaller. Based on these 
results, and based on the absence of significant disk thick- 
ening in our simulations, we are confident that Eq. (14) is 
sufficiently accurate for all realistic gas disks. 



APPENDIX E: VALIDITY CHECK OF THE 
REDUCED POISSON EQUATION FOR THE 
GAS DISK 

Throughout this paper we have assumed that the radial po- 
tential gradients of the disk are negligible compared to the 
vertical gradients, such that the Poisson equation reduces 
to Eq. (14). We now test this assumption by computing the 
ratio 

with $ g the gravitational potential of the gas disk. For a 
realistic, analytical disk model, our assumption will be valid 
as long as e « 1. 

Consider the Miyamoto & Nagai (1975) potential: 

$.(fl,s) = — — GMg = . (E2) 
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